C.................... RSLPSD.FOR ........................................... 
C--- subroutine loops is the main sequencing control program. it calls sub-
C--- progs DFIJ @ TFIJ to obtain the error functions FIJ for the density
C--- and temperature d.e.'s. it sets up the Jacobian matrix in subroutine
C--- JMATRIX and solves for the increments using BDSLV. increments from
C--- the solver are tested to ensure non -ve densities (modified steepest
C--- descent).
      SUBROUTINE LOOPS(S,NEQ,N,TI,ITAU,JBZB,SEC,MAXIT,ISOL
     >  ,JK,JCON,IWIND,RATIN,RATIS,HE)
      IMPLICIT DOUBLE PRECISION(A-H,L,N-Z)
      REAL SEC
      INTEGER ION,NEQ,NFLAG,LAXIT
      COMMON/SLV/DELTA(1806),RHS(1806),WORK(50000)
      DIMENSION N(4,401),TI(3,401),F(20),S(NEQ,1),DCRQ(2,401),V(2)
     > ,DUM(9),DCR(2),N1SAVE(401),HE(401)
      COMMON/VN/U(2,401),BG(401),BM(401),GR(2,401),R(2,401),SL(401)
      COMMON/ND/ON(401),HN(401),N2N(401),O2N(401),PHION(401),TN(401)
      COMMON/ALT/Z(401),DT,DH,THF,EPS,TF,ITF,JMAX,JMAX1,ITER,ION
      COMMON/SAV/NSAVE(2,401),TISAV(3,401),FY(2,401),UN(401),EHT(3,401)
      COMMON/LPS/SZA(401),EPSN,DC,IMOD,I1,I2,IDNDT,IPMX,IPP,ISKP
      DATA ICALL/3/,DCRT/0.0/,ZLBDY/90.0/,ZTBDY/90.0/
     > ,DCR/2*0.0/,ZLBSEC/0.0/

      JBT=JBZB     !... lower boundary point for Te and Ti
      JBN=JBZB     !... lower boundary point for O+ and H+ 

      !..... Get indices for the lower boundaries. Z(JBZB) is alt
      !..... from the specified absolute lower boundary
      IF(ZLBDY.LT.Z(JBZB)) ZLBDY=Z(JBZB)
      DO 41 J=1,JMAX/2
         IF(Z(JBN).LE.ZLBDY) JBN=J
         IF(Z(JBT).LE.ZLBDY) JBT=J
 41   CONTINUE

      !-- If F2 peak too close to equat point, switch off NmF2 normalize 
      IF((Z(JMAX/2+1)-Z(JK)).LT.50.0)  RATIN=-1.0
      IF((Z(JMAX/2+1)-Z(JCON)).LT.50.0)  RATIS=-1.0
      !-- Save current densities for dn/dt, dT/dt. FACIRI is a factor for 
      !-- scaling FLIP to the IRI nmF2. N1SAVE is used for restoring N(1,J)
      DO 43 J=1,JMAX
           !--- This section for scaling to IRI NmF2
           FACIRI=1.0
           ALPHA=5.0E-7
           IF(Z(J).GT.Z(JK)) ALPHA=1.0E-7
           IF(RATIN.GT.0.AND.J.LT.(JMAX+1)/2) 
     >         FACIRI=(RATIN-1.0)*EXP(-ALPHA*ABS(SL(J)-SL(JK)))+1
           IF(RATIS.GT.0.AND.J.GT.(JMAX+1)/2) 
     >          FACIRI=(RATIS-1.0)*EXP(-ALPHA*ABS(SL(J)-SL(JCON)))+1
           NSAVE(1,J)=1.0*(0.0+1.0*FACIRI)*N(1,J)

         N1SAVE(J)=N(1,J)
         NSAVE(2,J)=N(2,J)
         TISAV(1,J)=TI(1,J)
         TISAV(2,J)=TI(2,J)
         TISAV(3,J)=TI(3,J)
         DCRQ(1,J)=1.0
         DCRQ(2,J)=1.0
 43   CONTINUE

      !.. PREDIC was removed 19 June 2000: no longer makes sense with ExB

      !... start main loop here; if ifij=1 temps, if ifij=2 obtain densities
      DO 226 IFIJ=I1,I2
      IF(ISOL.EQ.0.AND.IFIJ.EQ.2) ITF=2
      IF(IFIJ.EQ.1) JB=JBT
      IF(IFIJ.EQ.2) JB=JBN
      JBS=JMAX-JB+1
      MIT=JMAX-2*JB+2
      IEQ=2*(MIT-2)
      MIT1=MIT-1
C
C.... avden is called to obtain average midpoint densities and temps
      IF(IFIJ.EQ.I2) CALL AVDEN(JMAX,TI,0,0.0D0,0.0D0)
C
C.... iter loop-- on each iteration the jacobian is formed and solved for
C.... the incs of density or temp. the density and temp are updated.
C.... if the densities at successive iterates are sufficiently close
C.... the solution has been found
C
      LAXIT=IABS(MAXIT)
      IF(DT.LE.120.AND.LAXIT.LT.29) LAXIT=29
      DO 220 ITER=1,LAXIT
C
C.... boundary conditions on density and temperature
      DO 321 J=1,JMAX
        !.. IF(Z(J).LT.ZTBDY) GO TO 321     2020-11-19
        IF(J.GT.JB.AND.J.LT.JBS) GO TO 321
          CALL HOEQ(1,J,JMAX,DT,THF,DUM,N,TI,NSAVE)
          N(1,J)=DUM(1)
          N(2,J)=DUM(2)
          IF(IFIJ.EQ.1) CALL HTEQCL(J,JMAX,N,TI)  !.. T boundary
 321   CONTINUE
C
C.... compute fij   profile and *s* matrix which bdslv uses to
C.... solve the linear system ((s))*(delta)= (fij)
C.... do loop for calling fij ....
C
       DO 145 J=2,MIT1
       KR=2*(J-2)
       JC=J+JB-1
       IF(IFIJ.EQ.1.AND.Z(JC).LT.0.000.AND.Z(JC).GT.200)
     >  WRITE(6,485) J,JC,INT(Z(JC)),TI(1,JC),TI(2,JC)
     > ,TI(3,JC)
 485  FORMAT(3I6,9F7.0)
       IF(IFIJ.EQ.1)  CALL TFIJ(JC,0,IPR,N,TI,F,JB,JBS,V,HE)
       IF(IFIJ.EQ.2)  CALL DFIJ(JC,0,IPR,N,TI,F,JB,JBS,V)
       RHS(KR+1)=F(1)
       RHS(KR+2)=F(2)
 145  CONTINUE
C147  CONTINUE
C
C.... ijm is a switch for calling jmatrix or not ....
       IJM=0
       IF((ITER/ICALL)*ICALL.EQ.ITER) IJM=1
       IF(ITER.LE.1) IJM=1
       IF(DCRT+DCR(1)+DCR(2).NE.3) IJM=1
C
       IF(IJM.EQ.1) CALL JMATRX(S,RHS,IEQ,IPR,DT,JMAX,N,TI,F,JB,JBS
     >   ,MIT,ITER,HE)
C
C
C.... invert the jacobian matrix *s* in the inversion routine *bdslv*.
C.... the increments are stored in array delta in this order
C.... x(1...n,j),x(1...n,j+1),x(1...n,j+2),....x(1...n,jmax-1)
C
      CALL BDSLV(IEQ,3,S,0,RHS,DELTA,WORK,NFLAG)
      IF(NFLAG.EQ.0)GO TO 55
      WRITE(6,215) ITER,NFLAG
      CALL RUN_ERROR    !.. print ERROR warning in output files
      STOP !3
C
C
 55     IDIV=0
        DCR(1)=1
        DCR(2)=1
        DCRP=1.0
        DCRT=1.0
        IF(IPR.NE.3) GO TO 135
C
C.... test ti increments 'dinc' to ensure ti>0 (mod steep descent) ......
       DO 137 I=1,IPR
       DO 137 J=2,MIT1
       DCRP=1.0
       ION=3
       IF(I.EQ.IPR) ION=2
       DINC=DELTA(2*J-ION)
C.... if dinc exceeds a fraction dc of ti set the factor dcrt so that ti>0
       JC=JB+J-1
       IF(ABS(DINC/TI(I,JC)).GT.EPSN) DCRP=DC*ABS(TI(I,JC)/DINC)
       IF(DCRP.LT.DCRT) DCRT=DCRP
C--- needs to be fixed for I>3       DCRQ(I,JC)=DCRP
 137   CONTINUE
 135   IF(IPR.NE.2) GO TO 144
C
C.... test n increments 'dinc' to ensure n>0 (mod steepest descent)
       DO  142 I=1,IPR
       DO 142 J=2,MIT1
       DCRP=1.0
       JC=JB+J-1
       DCRQ(I,JC)=1.0
       ION=3
       IF(I.EQ.IPR) ION=2
       DINC=DELTA(2*J-ION)
       IF(DINC.LE.0) GO TO 142
       IF(ABS(DINC/N(I,JC)).GT.EPSN) DCRP=DC*ABS(N(I,JC)/DINC)
       IF((DCRP.LT.DCR(I)).AND.(Z(JC).GT.0.0)) DCR(I)=DCRP
       IF(ITER.GT.0) DCRQ(I,JC)=DCRP
 142    CONTINUE
C
C
C... add iterative increment to the array 'n' or 'ti' and test for
C... convergence when idiv=0. separate dcrs are switched off
      ADCR=DCR(1)
      IF(DCR(2).LT.DCR(1)) ADCR=DCR(2)
 144   DO  42 I=1,IPR
       ION=3
       IF(I.EQ.IPR) ION=2
       DO 42 J=2,MIT1
C....... first temperatures ....
       JC=JB+J-1
       DINC=DELTA(2*J-ION)
       IF(IPR.EQ.2) GO TO 59
       TI(I,JC)=TI(I,JC)-DINC*DCRT
       IF(ABS(DINC/TI(I,JC)).GT.1E-3)   IDIV=IDIV+1
       GO TO 42
 59    CONTINUE
C........ now densities ..........
       N(I,JC)=N(I,JC)-DINC*ADCR
       IF(ABS(DINC/N(I,JC)).GT.1E-3)  IDIV=IDIV+1
 42    CONTINUE
C47    CONTINUE
C
C.... if dcr(i) is not close to 1 by the tenth iteration, give up.
       IF(DCR(1).LE.1E-3.AND.ITER.GT.10) GO TO 230
       IF(DCR(2).LE.1E-3.AND.ITER.GT.10) GO TO 230
C
C... test to see if convergence has occurred and write out Te heat diagnostics.
      IF(IDIV.EQ.0)  THEN
        !.. Convergence diagnostic for Temperatures. 2022-10-15
        IF(IFIJ.EQ.-999) THEN  ! set to IFIJ=1 to turn on print
          DO J=2,JMAX/2
            IF(Z(J).LT.200)  
     >     CALL TFIJ(J,4,IPR,N,TI,F,JB,JBS,V,HE)
          ENDDO
        ENDIF
        GO TO 224
      ENDIF
 220  CONTINUE
      GO TO 230
 224  CONTINUE

      !... Convergence in T or N 
      !... restore saved value for use in MINORA if it non-converges
      IF(IFIJ.EQ.I2)  THEN
         DO J=1,JMAX
           NSAVE(1,J)=N1SAVE(J)
         ENDDO

         !... blend in the region below the lower boundary. Also reset
         !... lower boundary altitude if decent time elapsed since problem
         IF(JBN.GT.JBZB)
     >     CALL SMOOTH(JMAX,JBN,JBZB,Z,N,NSAVE,DT,THF,DUM,TI)
         IF(SEC-ZLBSEC.GT.1800) ZLBDY=Z(JBZB)
      ENDIF
      IF(IFIJ.EQ.I2) RETURN
 226  CONTINUE

 
      !............ Deal with non-convergence
 230  CONTINUE
      IF(MAXIT.LT.0) GO TO 367
      IF(ITF.EQ.3.AND.IFIJ.EQ.2) GO TO 367
 
      !... restore old values for reduced time step
      DO 255 J=1,JMAX
        TI(1,J)=TISAV(1,J)
        TI(2,J)=TISAV(2,J)
        TI(3,J)=TISAV(3,J)
        N(1,J)=N1SAVE(J)
        N(2,J)=NSAVE(2,J)
 255  CONTINUE
C
      !.. Test to see if convergence difficulty below ZLBDY only
      IMOD=1         !... assume only below ZLBDY
      ABSBDY= 200    !... absolute boundary
      !... now test points above the boundary
      DO 257 J=JBN,JMAX+1-JBN
        IF(Z(J).LT.ABSBDY.AND.DCRQ(1,J).NE.1) IMOD=1
 257  CONTINUE
C
      IF(IMOD.EQ.0.OR.Z(JBN).GE.ABSBDY+1.0.OR.DT.GT.100) GO TO 365
      !.. If difficulty is only below ZLBDY km then raise lower boundary
      ZLBDY=(ABSBDY+2.0*Z(JBN))*0.333333
      ZLBSEC=SEC

      !... blend in low and high altitudes
      IF(IFIJ.EQ.2.AND.JBN.GT.JBT) 
     >   CALL SMOOTH(JMAX,JBN,JBZB,Z,N,NSAVE,DT,THF,DUM,TI)

      WRITE(6,373) NINT(ZLBDY)
 373  FORMAT(5X,'*** The O+ - H+ lower boundary has been raised,'
     > ,1X,'to',I4,' km ***')
C....... else return negative ITER to reduce time step
 365  ITER=-ITER
      IF(DT.GE.5) RETURN
 367  CONTINUE
C
C....  non-convergence diagnostics. grid points where problem occured
      IF(IFIJ.EQ.1) WRITE(3,115) ITER,DCRT,DCR(2),ZLBDY
      IF(IFIJ.EQ.1) WRITE(6,115) ITER,DCRT,DCR(2),ZLBDY
      IF(IFIJ.EQ.2) WRITE(3,116) ITER,DCR(1),DCR(2),ZLBDY
      IF(IFIJ.EQ.2) WRITE(6,116) ITER,DCR(1),DCR(2),ZLBDY
 115    FORMAT(/'  NON-CONVERGENCE IN TEMPERATURE SOLUTION, ITER='
     > ,I3,'  DCR1=',1P,E10.1,'  DCR2=',E10.1,' ZLBDY=',0P,F8.1)
 116    FORMAT(/'  NON-CONVERGENCE IN O+ - H+ SOLUTION, ITER='
     > ,I3,'  DCR1=',1P,E10.1,'  DCR2=',E10.1,' ZLBDY=',0P,F8.1)

      DO 535 J=JB+1,(JMAX+1)/2
      IF(DCRQ(1,J)+DCRQ(2,J).NE.2)
     > WRITE(6,537) J,INT(Z(J)),DCRQ(1,J),DCRQ(2,J)
 535  CONTINUE
 537  FORMAT(2X,'Northern hemisphere:- problem at point=',I5,' ALT=',I8
     > ,1P,4E8.1)
C
      DO 539 J=(JMAX+1)/2,JBS-1
      IF(DCRQ(1,J)+DCRQ(2,J).NE.2)
     > WRITE(6,541) J,INT(Z(J)),DCRQ(1,J),DCRQ(2,J)
 539  CONTINUE
 541  FORMAT(2X,'Southern hemisphere:- problem at point=',I5,
     >  ' ALT=',I8,1P,4E8.1)
      ITER=-1
      RETURN
 215   FORMAT('  ITERATION =',I5,'  RETURN FROM BDSLV =',I5)
C111    FORMAT(1H ,I4,1P,22E12.2)
C114    FORMAT(1X,1P,15E8.1)
      END
C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE HOEQ(ISW,J,JMAX,DT,THF,DUM,N,TI,NSAVE)
C*********  finds the new chemical equilibrium densities
C*********  of o+ and h+ at point j for fraction implicitness
      IMPLICIT DOUBLE PRECISION(A-H,L,N-Z)
      COMMON/ND/ON(401),HN(401),N2N(401),O2N(401),PHION(401),TN(401)
      DIMENSION RTS(199),N(4,401),TI(3,401),NSAVE(2,401),DUM(9)
       CALL RATS(J,TI(3,J),TI(1,J),TN(J),RTS)
C......... evaluate loss rates
      LOPLS=(RTS(3)*N2N(J)+RTS(4)*O2N(J))
      LOPLS2=+RTS(2)*HN(J)
      LHPLS=RTS(1)*ON(J)
      DUM(1)=PHION(J)/LOPLS + 1.0E-4   !.. stops O+ density going to zero
      DUM(2)=N(1,J)*LOPLS2/LHPLS
       RETURN
       END
C::::::::::::::::::::::::::: HTEQCL :::::::::::::::::::::::::::::::::::::
C.... lower boundary for Te using Newton to solve heating=cooling
      SUBROUTINE HTEQCL(J,JMAX,N,TI)
       IMPLICIT DOUBLE PRECISION(A-H,L,N-Z)
      COMMON/ND/ON(401),HN(401),N2N(401),O2N(401),PHION(401),TN(401)
      COMMON/VN/U(2,401),BG(401),BM(401),GR(2,401),R(2,401),SL(401)
      COMMON/SAV/NSAV(2,401),TISAV(3,401),FY(2,401),UN(401),EHT(3,401)
      DIMENSION N(4,401),TI(3,401),HLOSE(3),CLRT(3)
      JITER=0
      TI(1,J)=TN(J)
      TI(2,J)=TN(J)
      !... At night Te=Tn
      IF(EHT(3,J).LE.0.0) THEN
         TI(3,J)=TN(J)
         RETURN
      ENDIF
C..... Start a Newton procedure to solve local equilibrium      
C..... evaluate the heating and cooling rates for T and T+delT
      H=0.0001*TI(3,J)    ! delta T for Newton
 20   CONTINUE
        CALL GET_LOSS(2,J,TI,N,BM,HLOSE,CLRT)
        FEX=EHT(3,J)*1.6E-12/BM(J)-HLOSE(2)-CLRT(2)
        TI(3,J)=TI(3,J)+H
        CALL GET_LOSS(2,J,TI,N,BM,HLOSE,CLRT)
        DEX=(EHT(3,J)*1.6E-12/BM(J)-HLOSE(2)-CLRT(2)-FEX)/H
        TI(3,J)=TI(3,J)-H-FEX/DEX
C......... ensure te is physically reasonable
        IF(TI(3,J).LT.0.0) TI(3,J)=TN(J)
        IF(TI(3,J).GT.1.E+4) TI(3,J)=TN(J)
        JITER=JITER+1
        IF(JITER.GT.22) GO TO 30        !.. give up
      IF(ABS(FEX/(DEX*TI(3,J))).GT.1.0E-3) GO TO 20
      IF((J.NE.1).OR.(J.NE.JMAX)) RETURN
 30   CONTINUE
      IF(JITER.LT.22) RETURN
      !WRITE(6,*)'   JITER>22 IN HTEQCL J=',J
      TI(3,J)=TN(J)
        RETURN
        END
C::::::::::::::::::::::: JMATRX ::::::::::::::::::::::::::::::::::::
C......  JMATRX EVALUATES DELF/DELT USING STEPHANSONS METHOD.
C......  S(L,M) = DEL(FIJ) / DEL(N)
C......  WHERE N IS THE DENSITY OF ION IV AT ALTITUDE JV. L AND M
C......  ARE COMPUTED FROM JF...IV. THE ARRAY RHS CONTAINS VALUES OF
C...... FIJ SAVED FROM PREVIOUS COMPUTATION

      SUBROUTINE JMATRX(S,RHS,INEQ,IPR,DT,JMAX,N,TI,F,JB,JBS,MIT,ITER
     >  ,HE)
      IMPLICIT DOUBLE PRECISION(A-H,N-Z)
      DIMENSION RHS(1806),S(INEQ,8),N(4,401),TI(3,401),F(20),V(2)
     >  ,HE(401)
C
      DO 180 KZS=1,7
      DO 180 JZS=1,INEQ
 180  S(JZS,KZS)=0.0
C
      DO 80 JF=2,MIT-1
      J1=MAX0(2,JF-1)
      J2=MIN0(JF+1,MIT-1)
      DO 80 IV=1,2
      DO 80 JV=J1,J2
C
      L=2*(JF-2)
      IF(JF.LE.3)M=2*(JV-2)+IV
      IF(JF.GT.3)M=2*(JV-JF)+IV+4
      KRV=2*(JV-2)+IV
C
      JVC=JB+JV-1
      JFC=JB+JF-1
C
C... H FOR TI CALCULATIONS ..........
      IF(IPR.EQ.2) GO TO 68
      H=1.E-4*TI(IV+1,JVC)
      TI(IV+1,JVC)=TI(IV+1,JVC)+H
      CALL TFIJ(JFC,1,IPR,N,TI,F,JB,JBS,V,HE)
      TI(IV+1,JVC)=TI(IV+1,JVC)-H
       GO TO 70
 68    CONTINUE
C....... DENSITIES ........
      H=1.E-4*N(IV,JVC)
      N(IV,JVC)=N(IV,JVC)+H
      CALL DFIJ(JFC,1,IPR,N,TI,F,JB,JBS,V)
      N(IV,JVC)=N(IV,JVC)-H
 70     CONTINUE
      RH=1/H
      IF(JF.LE.3)  S(L+1,M)=(F(1)-RHS(L+1))*RH
      IF(JF.LE.3)  S(L+2,M)=(F(2)-RHS(L+2))*RH
      IF(JF.GT.3)  S(L+1,M-1)=(F(1)-RHS(L+1))*RH
      IF(JF.GT.3)  S(L+2,M-2)=(F(2)-RHS(L+2))*RH
C
   80 CONTINUE
      RETURN
      END
C:::::::::::::::::::::::::::: SMOOTH :::::::::::::::::::::::
C....... Smoothing out O+ profile when boundary above 120 km by 
C....... fitting an exponential from absolute lower boundary to the O+ 
C....... lower boundary. This profile is then blended with the
C....... local equilibrium values from subroutine HOEQ. The second
C....... call to HOEQ is to get the H+ density consistent with O+
      SUBROUTINE SMOOTH(JMAX,JBN,JBT,Z,N,NSAVE,DT,THF,DUM,TI)
      IMPLICIT DOUBLE PRECISION(A-H,L,N-Z)
      DIMENSION N(4,401),TI(3,401),NSAVE(2,401),DUM(9),Z(401)
        JINC=1
        JBS=JMAX-JBN+1
        JBTS=JMAX-JBT+1
        AL1=DLOG(N(1,JBN+JINC)/N(1,JBT))/(Z(JBN+JINC)-Z(JBT))
        AL2=DLOG(N(1,JBS-JINC)/N(1,JBTS))/(Z(JBS-JINC)-Z(JBTS))
C............ Northern hemisphere
        DO 511 J=JBT,JBN+JINC
           N(1,J)=N(1,JBT)*EXP(AL1*(Z(J)-Z(JBT)))
           CALL HOEQ(1,J,JMAX,DT,THF,DUM,N,TI,NSAVE)
           N(1,J)=(N(1,J)+DUM(1))*0.5
           CALL HOEQ(1,J,JMAX,DT,THF,DUM,N,TI,NSAVE)
           N(2,J)=DUM(2)
 511    CONTINUE   
        DO 512 J=JBT,JBN+10
C....           WRITE(73,417) J,Z(J),N(1,J),NSAVE(1,J)
 417       FORMAT(1X,I4,F6.1,1P,22E9.1)
 512    CONTINUE
C
C............ Southern hemisphere
        DO 515 J=JBS-JINC,JBTS
           N(1,J)=N(1,JBTS)*EXP(AL2*(Z(J)-Z(JBTS)))
           CALL HOEQ(1,J,JMAX,DT,THF,DUM,N,TI,NSAVE)
           N(1,J)=(N(1,J)+DUM(1))*0.5
           CALL HOEQ(1,J,JMAX,DT,THF,DUM,N,TI,NSAVE)
           N(2,J)=DUM(2)
 515    CONTINUE   
        DO 517 J=JBS-10,JBTS
C....           WRITE(74,417) J,Z(J),N(1,J),NSAVE(1,J)
 517    CONTINUE
        RETURN
        END
